#!/bin/bash
INPUT_MZML=$1
DATABASE=$2
DESIGN=$3
KIT_CORRECTION=$4

Help()
{
   # Display Help
   echo "Perform proteomic quantification of TMT16 sample via sage engine."
   echo
   echo "Syntax: sage_tmt16.sh [-h] -i INPUT_DIR -d DESIGN -k KIT_CORRECTION"
   echo "options:"
   echo "-h     Print this Help and exit"
   echo "-i     Directory path (no trailing /) containing all input mzML files"
   echo "-d     Openms-style design tsv"
   echo "-k     TMT16-kit correction matrix in ini file"
   echo
}

while getopts ":h" option; do
   case $option in
      h) # display Help
         Help
         exit;;
     \?) # incorrect option
         echo "Error: Invalid option"
         exit;;
   esac
done

echo "===== `date +%T` Running SageAdapter... (workflow 1/12) ====="

echo "INPUT_DIR = ${INPUT_MZML}"
echo "DATABASE = ${DATABASE}"
echo "DESIGN = ${DESIGN}"
echo "KIT_CORRECTION = ${KIT_CORRECTION}"

echo

export SAGE_LOG=trace

SageAdapter -in ${INPUT_MZML}/*/*.mzML -out out_0_sage.idXML \
    -threads 10 \
    -database "$DATABASE" \
    -min_len 6 \
    -max_len 40 \
    -min_matched_peaks 1 \
    -min_peaks 1 \
    -max_peaks 500 \
    -enzyme "Trypsin" \
    -precursor_tol_left -10.0 \
    -precursor_tol_right 10.0 \
    -fragment_tol_left -20.0 \
    -fragment_tol_right 20.0 \
    -variable_modifications 'TMTpro (K)' 'TMTpro (N-term)' 'Oxidation (M)' \
    -max_variable_mods 3 \
    -isotope_error_range 0,1 \
    -PeptideIndexing:IL_equivalent \
    -PeptideIndexing:unmatched_action warn \
    -debug 0 2>&1 | tee sage.log

echo "===== `date +%T` Running IDRipper... (workflow 2/12) ====="

mkdir -p 01sage
IDRipper -in out_0_sage.idXML -out 01sage -split_ident_runs
rm out_0_sage.idXML
for f in 01sage/*.idXML
    do
        mv "$f" "${f%.*}_sage.idXML"
    done

echo "===== `date +%T` Running Percolator... (workflow 3/12) ====="

mkdir -p 02percolator
for f in 01sage/*.idXML
do
OMP_NUM_THREADS=10 PercolatorAdapter \
        -in ${f} \
        -out 02percolator/`basename ${f} .idXML`_perc.idXML \
        -threads 10 \
        -subset_max_train 300000 \
        -decoy_pattern DECOY_ \
        -post_processing_tdc \
        -score_type pep \
        -debug 0
done

echo "===== `date +%T` Running IDScoreSwitcher... (workflow 4/12) ====="

mkdir -p 03IDScoreSwitcher
for f in 02percolator/*
do
IDScoreSwitcher \
        -in ${f} \
        -out 03IDScoreSwitcher/`basename ${f} .idXML`_fdr.idXML \
        -threads 2 \
        -new_score MS:1001491 \
        -new_score_orientation lower_better -old_score "Posterior Error Probability" -new_score_type q-value -debug 0
done

echo "===== `date +%T` Running IDFilter... (workflow 5/12) ====="

mkdir -p 04IDFilter
for f in 03IDScoreSwitcher/*
do
IDFilter \
        -in ${f} \
        -out 04IDFilter/`basename ${f} .idXML`_fil.idXML \
        -threads 2 \
        -score:pep "0.01" # the threshold here is 0.1 in quantms by default, change it to 0.01 during tweaking
done

echo "===== `date +%T` Running IsobaricAnalyzer... (workflow 6/12) ====="

mkdir -p 05isobaric
for f in ${INPUT_MZML}/*/*.mzML
do
IsobaricAnalyzer \
        -type tmt16plex \
        -in ${f} \
        -ini ${KIT_CORRECTION} \
        -threads 10 \
        -extraction:select_activation "auto" \
        -extraction:reporter_mass_shift 0.002 \
        -extraction:min_reporter_intensity 0.0 \
        -extraction:min_precursor_purity 0.0 \
        -extraction:precursor_isotope_deviation 10.0 \
        -tmt16plex:reference_channel 126 \
        -out 05isobaric/`basename ${f} .mzML`_iso.consensusXML \
        -debug 0
done 

echo "===== `date +%T` Running IDMapper... (workflow 7/12) ====="

mkdir -p 06IDMapper
for f in 05isobaric/*
do
IDMapper \
        -id 04IDFilter/`basename ${f} _iso.consensusXML`_sage_perc_fdr_fil.idXML  \
        -in ${f} \
        -threads 10 \
        -out 06IDMapper/`basename ${f} .consensusXML`_mapped.consensusXML \
        -debug 0
done

echo "===== `date +%T` Running FileMerger... (workflow 8/12) ====="

FileMerger \
        -in 06IDMapper/* \
        -in_type consensusXML \
        -annotate_file_origin \
        -append_method 'append_cols' \
        -threads 1 \
        -out ID_mapper_merge.consensusXML

echo "===== `date +%T` Running ProteinInference... (workflow 9/12) ====="

ProteinInference \
        -in ID_mapper_merge.consensusXML \
        -threads 10 \
        -picked_fdr true \
        -picked_decoy_string DECOY_ \
        -protein_fdr true \
        -Algorithm:use_shared_peptides true \
        -Algorithm:annotate_indistinguishable_groups true \
        -Algorithm:score_aggregation_method best \
        -Algorithm:min_peptides_per_protein 1 \
        -out ID_mapper_merge_epi.consensusXML \
        -debug 0

echo "===== `date +%T` Running IDFilter the 2nd... (workflow 10/12) ====="

IDFilter \
        -in ID_mapper_merge_epi.consensusXML \
        -out ID_mapper_merge_epi_filter.consensusXML \
        -threads 2 \
        -score:protgroup "0.01" -score:pep "0.01" -delete_unreferenced_peptide_hits -remove_decoys


echo "Running IDConflictResolver... (workflow 11/12)"

IDConflictResolver \
        -in ID_mapper_merge_epi_filter.consensusXML \
        -threads 5 \
        -out ID_mapper_merge_epi_filter_resconf.consensusXML

echo "===== `date +%T` Running ProteinQuantifier... (workflow 12/12) ====="

outname=`basename ${DESIGN} .tsv`

ProteinQuantifier \
        -method 'top' \
        -in ID_mapper_merge_epi_filter_resconf.consensusXML \
        -design ${DESIGN} \
        -out ${outname}_protein.csv \
        -mztab ${outname}_openms.mzTab \
        -peptide_out ${outname}_peptide.csv \
        -top:N 3 \
        -top:aggregate median \
        -top:include_all \
        -ratios \
        -threads 10 \
        -debug 0
